Numerical evolutions of nonlinear r-modes in neutron stars 
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Nonlinear evolution of the gravitational radiation (GR) driven instability in the r-modes of neu- 
tron stars is studied by full numerical 3D hydrodynamical simulations. The growth of the r-mode 
instability is found to be limited by the formation of shocks and breaking waves when the dimension- 
less amplitude of the mode grows to about three in value. This maximum mode amplitude is shown 
by numerical tests to be rather insensitive to the strength of the GR driving force. Upper limits 
on the strengths of possible nonlinear mode-mode coupling are inferred. Previously unpublished 
details of the numerical techniques used are presented, and the results of numerous calibration runs 
are discussed. 
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I. INTRODUCTION 

In recent years the gravitational radiation (GR) driven 
instability in the r-modes of rotating neutron stars has 
received considerable interest, both as a source of grav- 
itational waves for detectors such as the Laser Interfer- 
ometer Gravitational Wave Observatory (LIGO), and as 
an astrophysical process capable of limiting the rotation 
rates of neutron stars. In any rotating star, the r-modes 
are driven towards instability by GR ||: as the star 
emits gravity waves (primarily through a gravitomag- 
netic effect), the GR reaction acts back on the fluid by 
lowering the (already negative) angular momentum of the 
mode. This in turn causes the amplitude of the mode to 
grow. In most stars internal dissipation suppresses the 
r-mode instability, but this may not be the case for hot, 
rapidly rotating neutron stars For neutron stars 

with millisecond rotation periods, the timescale for the 
growth of the instability is about 40 s. In the absence of 
any limiting process, GR would force the dimensionless 
amplitude of the most unstable (to = 2) r-mode to grow 
to a value of order unity within about ten minutes of the 
birth of such a star. (At unit amplitude, the character- 
istic r-mode velocities are comparable to the rotational 
velocity of the star.) 

The strength of the GR emitted and the timescale 
on which the neutron star loses angular momentum and 
spins down depend critically on the maximum amplitude 
to which the r-mode grows. Initial estimates assumed 
that the amplitude would grow to a value of order unity 
before an undescribed nonlinear process saturated the 
mode. After saturation, it was assumed that the spin- 
down would proceed as a quasi-stationary process, re- 
ducing the angular velocity to one tenth of its initial 
value within about one year. In this scenario, gravita- 
tional waves from spindown events might be detectable 
with LIGO II g. 

However, at present no one knows with certainty how 
large the amplitude of the r-modes will grow. It may well 
be that the nonlinear hydrodynamics of the star might 
limit the growth of r-modes to very small values. This 



could happen, for instance, if the r-modes were to leak en- 
ergy by nonlinear coupling into other modes faster than 
GR reaction could restore it. In this case the r-mode 
instability would not play any interesting role in real as- 
trophysical systems. 

In a previous paper we presented the preliminary 
results of fully nonlinear, three-dimensional numerical 
simulations aimed at investigating the growth of r-modes. 
In our simulations, we modeled a young neutron star as 
a rapidly rotating, isentropic, Newtonian polytrope; we 
added a small-amplitude seed r-mode and we solved the 
hydrodynamic equations driven by an effective GR reac- 
tion force. We found that r-mode saturation intervenes 
at amplitudes far larger than expected (~ 3), supporting 
the astrophysical relevance of r-modes and the possibil- 
ity of detecting r-mode gravity waves. The details of 
the GR signature emitted by the r-mode instability that 
we observe in our simulations are rather different than 
previously envisioned, and these details suggest that this 
radiation may be more easily detected than previously 
thought: the radiation is more monochromatic and is 
emitted in a shorter, more powerful burst (see Ref. 
and the final Section of this paper). 

Our results are compatible with the conclusions of 
Stergioulas and Font , who performed relativistic simu- 
lations of r-modes on a fixed neutron-star geometry, and 
found no saturation even at large amplitudes. A second 
point of comparison can be made with the work of Schcnk 
and colleagues They have attacked the problem an- 
alytically, developing a perturbative formalism to study 
the nonlinear interactions of the modes of rotating stars, 
and proving that the couplings of r-modes to many other 
rotational modes are small (they are forbidden by selec- 
tion rules, or they vanish to zeroth order in the angular 
velocity of the star). 

The present work is meant as a complement to Ref. : 
throughout these pages, we describe our simulations in 
greater detail; we discuss their relevance and their limita- 
tions in the light of Refs. @, ||; and we present the results 
of several additional simulations aimed at enlightening 
particular aspects of the problem. In Sec. II we write 
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down the basic hydrodynamic equations, and we define 
a number of mathematical quantities that will be used 
to monitor the nonlinear evolution of the r-modes. In 



Sec. [I] we implement the effective current-quadrupole 
gravitational radiation-reaction force. In Sec. IV we 
integrate the fluid equations with r-mode initial data 
in slowly rotating stars, and we compare the results 
with the small-amplitude, slow-rotation analytical ex- 
pressions: we demonstrate that the integration repro- 
duces faithfully the analytical p redic tions to the expected 
degree of accuracy. In Sees. ^VIII we study the nonhn- 
ear evolution of r-mode initial data in rapidly rotating 
stars, concentrating on the nonlinear saturation of the 
r-modes, and analyzing in detail the evolution of several 
hydrodynamical quantities. Finally, we summarize our 



conclusions in Sec. IX 



II. BASIC HYDRODYNAMICS 

We study the solutions to the Newtonian fluid equa- 
tions. 



dtp + V- (pv) = 0, 



(1) 



p{dtv + v- Vv)^-Vp- pV$ + pF^^, (2) 



dtT + V ■ {tv) = 0, 



(3) 



where v is the fluid velocity, p and p are the density and 
pressure, $ is the Newtonian gravitational potential, and 
F'^^ is the gravitational radiation reaction force. Equa- 
tion (||) is a recasting of the energy equation for adiabatic 
flows, where r is the entropy tracer for polytropic 
equations of state, r is related to the internal energy 
(per unit mass) e by the relation t = (ep)^/^, where 7 
is the adiabatic exponent. The Newtonian gravitational 
potential is determined by Poisson's equation, 



= AnGp, 



(4) 



while the gravitational radiation reaction force will be 
discussed in Sec. III. 

We solve Eqs. numerically in a rotating refer- 

ence frame, using the computational algorithm developed 
at LSU to study a variety of astrophysical hydrodynamic 
problems The code performs an explicit time in- 

tegration of the equations using a finite-difference tech- 
nique that is accurate to second order both in space and 
time, and uses techniques very similar to those of the fa- 
miliar ZEUS code [|ll]. For most of our simulations, we 
adopt a cylindrical grid with 64 cells in the radial direc- 
tion, and 128 cells in the axial and azimuthal directions. 

In the limit of slow rotation, we define the r-modes 
of rotating Newtonian stars (using the normalization of 
Lindblom, Owen and Morsink |3| ) as the solutions of the 



perturbed fluid equations having the Eulerian velocity 
perturbation 



(5) 



where R and f2o are the radius and angular velocity of 
the unperturbed star, ao is the dimensionless r-mode am- 
plitude, and Yif is a vector spherical harmonic of the 
magnetic type, defined by 

Y^^^W + l)]-'/'rVxirVYi^). (6) 

The r-mode frequency is given by 

(/-!)(/ + 2) 



l + l 



-fin 



(7) 



To monitor the nonlinear evolution of the r-modes, it 
is helpful to introduce nonlinear generalizations of the 
amplitude and frequency of the mode. These quantities 
are defined most conveniently in terms of the current 
multipole moments of the fluid, 



Ji^= I pr'v-YiZ*Sx. 



(8) 



In slowly rotating stars, the J22 moment is proportional 
to the amplitude of the m = 2 r-mode, the most unsta- 
ble mode, and the one that we will study. To track the 
evolution of this mode even in the nonlinear regime, we 
define the normalized, dimensionless amplitude 



2IJ22I 



(9) 



where M is the total mass of the star and J is defined by 
1 



in 



pr^d^x ~ / prVr. (10) 



The quantity J is evaluated once and for all at the begin- 
ning of each of our evolutions. For slowly rotating stars, 
the definition (^) of the mode amplitude reduces to the 
one given by Eq. (||). 

In slowly rotating stars, and in all situations where 
the leading contribution to J22 comes from the m ~ 2 r- 
mode, the time derivative dJ22/dt is proportional to the 
frequency of the mode: dJ22/dt = iujJ22- Thus we are 
led to define the nonlinear generalization of the r-mode 
frequency as 



1 



\J2 



dJ: 



22 



dt 



(11) 



As shown by RezzoUa, et al. ||lj], we can re-express 
dJ22/dt as an integral over the standard fluid variables. 



(1) _ 



22 



dt 



v-{VY2T)-v-y<^-Y2T d'x. (12) 
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The definitions, Eqs. (|9|) and (|11|), of mode amphtude 
and mode frequency are very stable numerically, because 
they are expressed in terms of integrals over the fluid 
variables. In the Appendix, we give explicit expressions 
for J22 and j'^ in the cylindrical coordinate system used 
in our numerical analysis. 

While we monitor the nonlinear evolution of the r- 
mode, we are also interested in tracking the star's average 
angular velocity as well as its degree of differential rota- 
tion. With this in mind, we define the average angular 
velocity 



n EE J//, 



(13) 



where the angular momentum and the moment of inertia 
are given respectively by 



J = y pn7'^n{n7, z,ip) (fix, 



2 j3 
pzu a X 



(14) 



(15) 



Here w is the cylindrical radial coordinate, and the lo- 
cal angular velocity ^[vj, 2, ip) = v^/vj, where is the 
proper azimuthal component of the fiuid velocity. We 
also define the average differential rotation Ail as the 
weighted variance of fi, 



I pzu^(n-nYd^x 



(16) 



III. RADIATION-REACTION FORCE 

The gravitational radiation-reaction force due to a 
time-varying current quadrupole is given by the expres- 
sion 



-Vj<^aUXkS\f - eaMXkXmS\^ , (17) 

see Blanchet and Eq. (20) of RezzoUa et al. @. 

Here Sj^^ represents the n'th time derivative of the cur- 
rent quadrupole tensor, 

Sjk = j p{xx v)(jXk)(fix\ (18) 

tjki is the totally antisymmetric tensor, and the vector Xk 
represents the Cartesian coordinates of the point at which 
the force is evaluated. The parameter k that appears 
in Eq. (^7|) has the value k = 1 in general relativity. 
For reasons discussed below, we find it useful to consider 
other values of k, in our numerical simulations. 



We find that a straightforward application of Eq. (|lj^ 
in numerical evolutions is nearly impossible. There are 
two problems: first, it is very hard to evaluate reliably 
time derivatives of such a high order; second, various 
sources of numerical noise (even small errors in the initial 
equilibrium configuration of the fiuid, and the numerical 
drift of the center of mass) can generate contributions to 
the current quadrupole tensor that overwhelm those of 
the pure r-mode motion. So we need to introduce spe- 
cial numerical techniques and simplifications to overcome 
these problems. 

In order to reduce the infiuence of extraneous noise 
sources on the evolution, it is helpful to re-express the 
current quadrupole tensor in terms of the current multi- 
pole moments defined in Eq. (||). There is a one-to-one 
correspondence between the J2m current multipoles and 



^yy ^xx ~t- 2%Sxy — 




(19) 
(20) 
(21) 



In a slowly rotating star, the ni ~ 2 r-mode excites J22, 
but not J21 and J20. In contrast, the principal sources 
of numerical noise contribute primarily to J20. Thus, 
we evaluate only the J22 contribution to F^^: we use 
Eq. ([l7|) to evaluate taking the Sij determined from 
Eqs.^l^)-(|l|), but setting J21 = J20 = 0. We find that 
this scheme reduces considerably the numerical noise in 
the radiation reaction force, and reproduces faithfully the 
analytic al de scription of r-modes in slowly rotating stars 
(see Sec.0. 

The second major problem is evaluating the numerical 
time derivatives of Sjk, or equivalently the time deriva- 
tives of J22 . Whenever the radiation-reaction timescale is 
much longer than the r-mode period 27r/a;, the dominant 

(n) 

contribution to the derivatives S)J comes from terms 
proportional to powers of the r-mode frequency: 



c(") 



jk, 



7(") 



22- 



(22) 



Even when the r-mode amplitude becomes large, the ex- 
pression (|2^) will be accurate as long as the timescale 
for the evolution of a and lu is longer than 27r/a;. Now, 

J22 and J22' are easily evaluated using the integral ex- 
pressions in Eqs. ^ and (p^; thus, the time derivatives 
needed in Eq. (17) are given simply by J22'' = ^^'^22' 



and 722' = —uj'^J22, where we determine lu numerically 
using Eq. (|l^). In the Appendix, we present explicit 
expressions for the components of the effective radiation- 
reaction force in cylindrical coordinates. 
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TABLE I: Physical parameters for the equilibrium models. 



Parameter 


Symbol 


Slow 


Fast 






CI, C2 


C3-C8 


polytropic index 


n 


1 


1 


total mass 


M 


1.4 M0 


1.4 M0 


equatorial radius 




12.7 km 


18.4 km 


polar R/equatorial R 




0.98 


0.59 


nonrotating" R 


-Ro 


12.5 km 


12.5 km 


angular velocity 


no 


1.45 krad/s 


5.34 krad/s 


rotation period 


Po 


4.32 ms 


1.18 ms 


energy ratio 


T,ot/\W\ 


3.98 X 10"^ 


0.10072 


simulation RR timescale 


'rr 


0.459 Po 


9.4374 Po 


physical RR timescale 


'rr. 


2.8 X lO'^Po 


4.2 X 10*Po 



"i?0 is the radius of the nonrotatine star of the same mass. For 
n = 1 polytropes Rq = nK/2G |ll6|| where K is the polytropic 
constant. 



0.16 ■ 




0.05 0.10 0.15 0.20 
t/Po 

FIG. 1: Evolution of the r-mode amplitude in a slowly ro- 
tating star. The solid curves plot the results of numerical 
evolutions (with and without gravitational radiation reac- 
tion) while the dashed curves plot the analytical predictions. 
For the curves marked "free," k = 0; for the curves marked 
"forced," k ~ 6 x 10^ 



IV. CALIBRATION RUNS 

In order to test the accuracy of our hydrodynamic evo- 
lution code and of our approximations for the gravita- 
tional radiation-reaction force, we investigate the evolu- 
tion of a small-amplitude r-mode in a slowly rotating 
star. 

We provide initial data for this study by solving the 
time-independent fluid equations for a slowly, rigidly ro- 
tating stellar model. We model the neutron star as an 
n = 1 polytrope, generated by the self consistent field 
technique developed by Hachisu . Table | shows the 
physical parameters for this model, labeled Slow; in par- 
ticular, the ratio of rotational kinetic energy to gravi- 
tational binding energy is Trot/|W^| = 0.00398, and the 
angular velocity is 26% of the maximum possible value 
(estimated as rimax = ^y/irGp). 

We then adjust the velocity field of this equilibrium 
model by adding the velocity perturbation of an m = 2 
r-mode of amplitude ag : 

v = mnne^ + anmo(^^yRe{Yi). (23) 

In the Appendix, we write out explicitly the components 
of this initial velocity field in our cylindrical coordinate 
system. Because Eq. ( p3|) is the exact representation 
of a pure m = 2 r-mode only in the small-amplitude, 
small- rotation limit, we expect that the frequency and 
the amplitude measured in our numerical experiment us- 
ing Eqs. (|9|)-(|ri|) will be different from their theoreti- 
cal values Q!o and — |f^o by terms of order 0(q:^), and 

We perform two numerical integrations of the equa- 
tions of motion, for this slowly rotating initial configura- 
tion. In the first run (CI), we let the star evolve under 
purely Newtonian hydrodynamics, setting the strength k 
of the radiation reaction (RTh to zero. In the second run 



(C2), we force the mode by setting k ~ 6 x 10^. With 
this unphysically large value the amplitude of the r-mode 
grows appreciably within a time that we can conveniently 
follow numerically. (The Courant limit for the evolution 
timestep is set by the speed of sound in the fluid, and by 
the size of the grid cells; for fig = 0.26riinax, one com- 
plete rotation of the star takes about 70000 timesteps). 

Figure |l| illustrates the evolution of the mode ampli- 
tude a in runs CI and C2, as a function of I/Pq, where 
Pq — 27r/rio is the initial rotation period of the star. 
The solid curves trace the numerical evolution of a [as 
defined in Eq. (^], whereas the dashed curves trace the 
theoretical predictions for this evolution, obtained in the 
small-amplitude, slow-rotation limit 

When K = 0, the theoretical prediction for the evolu- 
tion of the amplitude is just a = ao, and we can verify 
that the numerical evolution tracks the analytical curve 
within the expected deviations of order a^. When k ^ 0, 
the analytical prediction for the evolution of the ampli- 
tude is 

a^aoe*/^"^", (24) 
where the radiation-reaction timescale is given by 

— = '^''(^X Aj^'I^'^'o- (25) 
TGR V405y c7 

For this model, tgr = 0.46Po- As we can see in Fig. |^, 
the numerical evolution tracks the small-amplitude, slow- 
rotation analytical result within the expected accuracy, 
even if the radiation-reaction force is so unphysically 
strong. 

Although this slow rotation numerical evolution was 
only carried out over a small fraction (0.2) of a rotation 
period (and therefore over a small fraction of the r-mode 
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0.05 0.10 0.15 0.20 
tIPQ 

FIG. 2: Real and imaginary parts of the current multipole mo- 
ment J22 (in arbitrary units), for a slowly rotating star evolved 
without gravitational radiation reaction (run CI). The solid 
curves trace the numerical evolutions, the dashed curves trace 
the analytical predictions. 
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o 
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® -1.330 ■ 



-1.335 \ ^ ^ 

0.05 0.10 0.15 0.20 
tlPo 

FIG. 3: Frequency of the m = 2 r-mode, for a slowly rotat- 
ing star evolved without gravitational radiation reaction (run 
CI). The solid curve is determined numerically from Eq. ([l 
For comparison, the dashed line shows the analytical value 
ujo = — |f2o. 



oscillation period), the evolution extended for about 7.3 
dynamical times and 4.6 sound-crossing times. 

In Figs. ^ and ||, we display two additional diagnostics 
for the undriven (k — 0) slow-rotation evolution (CI). 
In Fig. ^ we plot the real and imaginary parts of the 
current multipole moment J22 : the solid curves trace the 
numerical evolution, whereas the dashed curves trace the 
analytical expression 



J22 ^ ^aMR^jnae^'^K 



(26) 



Again the deviations are within the expected accuracy 
of the analytical results. The deviations appear to be 
caused by the excitation of modes other than the pure 
m = 2 r-mode; the spurious excitations appear because 
the initial data [Eq. (^3|)] are only accurate to first order 
in a. Figure ^depicts the evolution of the frequency ui [as 
defined in Eq. (Kw]. The deviations from the analytical 
result, ivq — — ^^^o; are within the expected accuracy. 
The magnified scale used to display uj in Fig. |^ makes the 
presence of the small- amplitude, short-period extraneous 
modes quite apparent. 



T,ot/\W\ = 0.10072, and the angular velocity is 95% of 
its maximum value. 

We perform a numerical integration of the equations 
of motion starting from the rapidly rotating initial con- 
figuration. Fast, using Eq. (^) to add a slow-rotation, 
small-amplitude r-mode field, with ao — 0.1. Because 
the radiation-reaction force is so much stronger for this 
model (it is proportional to oc f2^), we find that we 
can set k = 4487, which yields an r-mode growth time 

(s) 

Tj^P — 9.43P0. This choice of k is still much larger than 
its physical value (unity) , but it should yield a reasonable 
picture of the nonlinear evolution of the r-mode, if the 
timescales for all the relevant hydrodynamical processes 
(including nonlinear couplings to other modes) are com- 
parable to Pq, or shorter. Indeed, if the average sound- 
crossing time Ts is representative of the relevant hydro- 
dynamical timescales, then our condition is satisfied: a 
rough estimate gives rg = Ro/cs — O.I6P0 — '''rr/60, 
where we have approximated eg as the average speed of 
sound in the equivalent spherical polytrope. 



A. Evolution of the r-mode amplitude 



V. EVALUATING THE SATURATION 
AMPLITUDE 

In our production runs, we investigate the nonlin- 
ear behavior of the r-mode in a rapidly rotating stellar 
model, under a variety of physical conditions (different 
initial amplitudes, and different values for the radiation- 
reaction coefficient k). 

Again, we provide initial data by solving the time- 
independent fluid equations for an n = 1 polytrope. 
The physical parameters for this model, labeled Fast, 
are reported in Table in particular, the ratio of ro- 
tational kinetic energy to gravitational binding energy is 



We follow the evolution through t = SSPq. Because the 
rotation of the star is progressively reduced by radiation 
reaction, and because the star develops differential rota- 
tion, at the end of the evolution the star has performed, 
on the average, only about 31 rotations [we obtain this 
number from J Cldt/ (27r)]. In Fig. ^ we plot the numeri- 
cally determined evolution of Re J22 ■ the curve is a very 
smooth sinusoid, whose frequency is essentially constant, 
and whose envelope is determined by the (relatively slow) 
evolution of the r-mode amplitude. So the approxima- 



tions used to compute u and J22^ (discussed in Sec. |I|) 
are in fact quite good in this situation. 

In Fig. ^ we plot the numerical evolution of the r- 
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FIG. 4: Numerical evolution of Re J22 (arbitrary units) for 
a rapidly rotating star driven by gravitational radiation re- 
action (production run C3). The sinusoidal approximation, 
used to compute uj and Jjj ' , is evidently appropriate for this 
run. 




4 8 12 16 20 24 28 32 

t/Po 

FIG. 5: Numerical evolution of the r-mode amplitude a in 
production run C3. 



mode amplitude a. At the beginning of the evolution, the 
computed diagnostic a agrees with the theoretical value 
ao to ^ 10%, within the expected accuracy. The growth 
is exponential (as predicted by perturbation theory) until 
a « 1.8. Then some nonlinear process begins to limit the 
growth, until the amplitude peaks at a = 3.35 and then 
falls rapidly within a few rotation periods. After this the 
r-mode is effectively not excited. 



B. A mechanism for r-mode saturation 

What nonlinear process is responsible for the behav- 
ior of the r-mode amplitude? What causes the mode to 
saturate and disappear from the star? To answer these 
questions, we study the evolution of the total mass, total 
angular momentum, and total kinetic energy of the star, 
which are plotted in Fig. ||. 

Because the mass is constant, the damping of the r- 



4 8 12 16 20 24 28 32 

t/Po 

FIG. 6: Evolution of total mass, total angular momentum, 
and total kinetic energy in production run C3. The quantities 
are plotted as fractions of their initial value. 



mode cannot be caused by ejection of matter from the 
simulation grid. On the other hand, we expect that the 
star should lose energy and angular momentum as it radi- 
ates gravitational radiation in accord with the prediction 
of general relativity llj, : 



H 
2 



1287r G 
225 1? 



K^^|J22P. (27) 



The evolution of the angular momentum mirrors this 
equation quite closely (within a few percent); for energy, 
however, Eq. ( p7| ) is only accurate until shortly after the 
catastrophic fall of the r-mode amplitude (at t ~ 28Po; 
see Fig. ^. Before that time, the star loses about 40% 
of its initial angular momentum and 36% of its initial ki- 
netic energy. After that time, the amplitude and (there- 
fore) the radiation-reaction force are much reduced, so 
J becomes essentially constant; however, the kinetic en- 
ergy continues to decrease, losing an additional 12% of 
its initial value during the next three rotation periods. 

If the r-mode were damped by a hydrodynamical pro- 
cess that conserved energy, such as the transfer of energy 
to other modes, then Eq. ( |27| ) should portray accurately 
the evolution of the kinetic energy. But this is not what 
we see: instead, some purely hydrodynamic process con- 
tinues to decrease the energy (by a sizable amount!) after 
the gravitational-radiation losses become negligible. 

We believe that we have identified this process. To 
first order in the amplitude, the r-mode is only a velocity 
mode; to second order, however, there is also an asso- 
ciated density perturbation, proportional to ¥32, which 
appears as a wave with four crests (two in each hemi- 
sphere) on the surface of the star. (We will present a 
quantitative analysis later in this section.) As the am- 
plitude reaches its maximum, these crests become large, 
breaking waves: the edges of the waves develop strong 
shocks that dump kinetic energy into thermal energy. In 
doing so they damp the r-mode. Figure ^ illustrates the 
surface waves at t = 28Po and t = 29Po along selected 
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4 8 12 16 20 24 28 32 
tIPQ 

FIG. 7: Theoretical and numerical evolution of the total en- 
ergy for production run C3. The total energy is plotted in 
units of the initial rotational energy (because the system is 
bound, Etot must be negative). 




FIG. 8: Isodensity surfaces showing breaking waves near the 
end of production run C3. 



meridional slices. 

Our code is written in such a way that the evolution 
of the shocks is always kinematically correct. In partic- 
ular, the continuity equation ensures that mass is prop- 
erly conserved, and that the proper density jump occurs 
across the shocks; and the Euler equation ensures that 
momentum is properly conserved, and that the proper 
pre- and post-shock velocities arise. However, our code 
does not allow the energy dissipated in the shock to in- 
crease the entropy of the fluid, which remains always 
barotropic and isentropic. Consequently, the presence of 
shocks shows up as a decrease in the total energy of the 
star. Indeed, in this production run (C3) gravitational 
radiation reduces the total energy by 9% of its initial 
magnitude through t — 28Pq, and dissipation in shocks 
subtracts a further 3% in the last three rotations. (For 
comparison, 12% of the initial kinetic energy is lost in the 
last three rotations compared to 36% before that time.) 

Ignoring the thermal effects of shocks is useful to re- 
duce the computational burden and the complexity of the 




0.2 0.4 0.6 0.8 1.0 
r/R 

FIG. 9: Radial amplitude density a{r) of the m = 2 r-mode 
for production run C3 at t = 22.5Po. 



hydrodynamic code, and it is in fact a fairly reasonable 
approximation for neutron star matter, where the pres- 
sure comes mostly from the Fermi pressure of the degen- 
erate neutrons, so the equation of state can be effectively 
modeled as temperature-independent. 

C. Radial structure of the r-mode amplitude 

We define the radial amplitude density a{r) (where r 
is the spherical radius) by expressing the integral Eq. (|^) 
for J22 in spherical coordinates, and omitting the radial 
integration: 

Q,(r)e«'>W = — / pr'^v ■Yjl*r'^ sin ededif. (28) 

JMB?VLq J 

We removed the absolute value around the integral 
for J22 so that we can keep track of the local 
mode phase (/'(r). With this definition, Q;exp[i0] = 
/ a{r) exp[i(/)(r)]dr/i?, where is the global phase of the 
r-mode. 

The amplitude air) is plotted in Fig. |9|for the produc- 
tion run C3 at the time t — 22.5Po- Throughout the en- 
tire evolution, the mode is concentrated mostly between 
the spherical radii r = 0.5i? and 0.9i?, and the shape of 
a(r) is fitted reasonably well by taking 6v cx {r/R)'^ [see 
Eq. (^)] and p cx {simrr / R) / (nr / R) as appropriate for a 
spherical, n = I polytrope. 

It is also interesting to study the phase coherence of 
the r-mode, which we define as 

, ra(r)|e*'^W -e^'^l^dr , ^ 

(A0)2 = J ^ ^1 \—. (29) 

J Q;(r)dr 

Figure |l^ plots the evolution of A^, which is small un- 
til the r-mode saturates at t ~ 26Po- For A(j) « 1, 
the local phase (/)(r) spans approximately 27r: the mode 
has lost coherence completely. In this situation, there 
are large regions in the star where the radiation-reaction 
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FIG. 10: Evolution of the phase-coherence function A<^ in 
production run C3. 
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FIG. 12: Numerical evolution of the average stellar angular 
velocity Q, in production run C3. 
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FIG. 11: Numerical evolution of the r-mode frequency uj in 
production run C3. 



force pushes out of phase with the local mode oscillations; 
this mismatch accelerates the damping of the mode. 



mentum decreases by 40%. (As the star spins down, it 
becomes less flattened, and the change in the moment of 
inertia accounts for the difference between the decrease of 
J and that of Cl.) The stability of the r-mode frequency 
has important implications for the possible detection of 
r-mode gravity waves (see Sec. p^ ). 

We also point out that the approximate expressions for 
the GR reaction force, Eqs. ([l7|)-(^, that we use here 
are accurate only when the motion of the fluid has nearly 
sinusoidal time dependence. Figure |ll| illustrates that 
the evolution in our simulation remains quite sinusoidal 
until about t — 28Po- After this point our expression for 
the GR reaction force is not reliable. After this point in 
our simulation, however, the fluid evolution is dominated 
by nonlinear hydrodynamic forces including shocks, and 
the GR reaction force is negligible. Thus our inability to 
model accurately the GR force during the late stages of 
the evolution does not effect our results. 



E. Growth of differential rotation 



D. Evolution of the r-mode frequency 

Figure |ll| shows the numerical evolution of the r-mode 
frequency w. The evolution of lo is quite smooth when 
the amplitude of the r-mode is large; when the ampli- 
tude is small (for t ^ lOPg and for t > 28Po), we see that 
other modes make noticeable contributions to J22, and 
therefore to w. At the beginning of the run, the numer- 
ical u) matches the theoretical prediction to within the 
expected accuracy of about 10%. These values for the 
frequency are also consistent with those obtained via a 
Fourier transform of J22 [@. 

Surprisingly, the r-mode frequency remains approxi- 
mately constant throughout the evolution, and it does 
not follow the decline of the average angular velocity 
(plotted in Fig. |l^). Altogether, the angular velocity 
decreases by about 22.5% while the total angular mo- 



During this simulation (run C3), the average differen- 
tial rotation Afi [defined in Eq. (pj|)] grows to a max- 
imum of approximately 0.4lf2 (see Fig. |l^). After a 
rapid increase in the first three rotation periods, when the 
linear r-mode eigenfunction of Eq. ( p3| ) evolves into its 
proper nonlinear, rapid-rotation form, AfJ/fJ increases 
approximately as a" ''^ until a ~ 1, and then approxi- 
mately as a until a begins to saturate. When a is maxi- 
mum. Ail = 0.25rj. As the amplitude falls, AO continues 
to grow (even more steeply) , as long as there is significant 
gravitational radiation. After t = 28Poi decreases to 
about 80% of its peak value. So the final configuration of 
the star (where the presence of the r-mode is essentially 
negligible) still has a very large differential rotation. 

But we should not concentrate exclusively on the av- 
eraged quantity AJl/n, which does not capture fully the 
spatial structure of differential rotation. Figure pi illus- 
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FIG. 13: Numerical evolution of the differential rotation Af2 
in production run C3. 




FIG. 14: Meridional structure of the differential rotation 
in production run C3. The plot shows the value of the 
azimuthally averaged angular velocity ^l(vj,z)/Cl, at time 
t = 25.6Po. 



trates the spatial dependence of the azimuthally averaged 
angular velocity, 

i7(ti7, = — j n{ru, z,(p) d(p, (30) 

at the time when the amplitude is maximum, t = 25.6Po- 
The differential rotation is confined mostly to a thin shell 
of material near the surface of the star, and is particu- 
larly concentrated near each polar cap. The bulk of the 
material in the star remains fairly rigidly rotating. 



F. Consistency of the radiation-reaction force 



followed by mass octupole and current quadrupole. To 
verify that our approximation is justified for the phys- 
ical states considered here, we evaluate the additional 
energy that would have been lost to gravitational waves 
throughout our simulation if we had included the lowest- 
order mass multipole terms. 

The mass multipole moments Qim are defined by 



Qhn = pr Yi*^d 



(31) 



In the presence of density oscillations with sinusoidal 
dependences in the coordinates t and Lp {i.e., Spim oc 
e'"''"*^'""^) the flux of energy into gravitational waves is 
given by |9[ 



8tt G 
5^ 

iTT G 



75 c° 



6615 c 



(32) 
(33) 



where Q2m and Qzm are, respectively, the mass 
quadrupole and mass octupole moments induced by these 
density fluctuations. Contributions of higher order are 
suppressed by very small fractional coefficients. 

Comparing Eqs. ( ^ ) and (^) with Eq. ( |27| ) we find 
that the contribution of the density oscillations associ- 
ated with the r-mode at frequency lo to the energy flux 
is negligible whenever 



3c2 IQ 



2m I 



16 J' 



« 1, 



5c^2 IQ 



3m I 



2m I 



2352 \J: 



< 1. 



2m I 



(34) 



We find that in our simulation both ratios are of or- 
der 10"'^ before the r-mode saturates (at t ~ 25Po)- 
The strongest contribution to the quadrupole term comes 
from Q22) although the Fourier transform of this moment 
does not show any definite frequency of oscillation. The 
strongest contribution to the octupole term comes from 
the Yi,2 dependence of the density in the m = 2 r-mode 
(see the next subsection). 

Between t ~ 25Po and t ~ 32Po (when a is back to 
its initial value ~ 0.1) the mass quadrupole term would 
have provided a correction of order 10% to the current 
quadrupole; although even then we see no evidence of a 
definite oscillation frequency correlated to the r-mode. 
Only after t ~ 32Po, when the fluid motion in the star 
becomes quite turbulent and the r-mode is very weak, is 
the gravitational radiation generated by the mass multi- 
poles comparable to the radiation from J22- 

On the whole, we find that our approximation which ig- 
nores the contributions from the mass multipoles is well- 
justified throughout the more interesting part of the evo- 
lution. 



In these simulations we have assumed that the only rel- 
evant contribution to the radiation reaction force comes 
from the current quadrupole moment, and in particular 
from J22. However, in the post-Newtonian approxima- 
tion to general relativity, the lowest-order contribution to 
radiation reaction comes from the mass quadrupole term. 



G. Density oscillations and mode saturation 

The evolution of the isodensity surfaces in our neu- 
tron star shows very clearly the presence of the lowest- 
order Eulerian density perturbation Sp associated with 
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FIG. 15: Numerical evolution of the mass moment (^32 (solid 
line) and of the r-mode amplitude a, in production run C3. 
The curve for Q32 was renormalized to emphasize the linear 
relation between a and Q32 during the growth of the r-mode. 



the m = 2 r-mode. The lowest order expression for 5p 
was derived by Lindblom, Owen and Morsink in the 
small-amplitude, slow-rotation approximation. Solving 
Eq. (5) of Ref. [g| with m = 2 and with polytropic index 
n — 1, and then substituting 6'^ back into Eq. (4) of 
Ref. I), we get 

where ja is the spherical Bessel function. The mass mul- 
tipole associated with this 5p is 



5Q 



32 



"°T5V3" 



7tt /2 nlRl 



(36) 



where jii?:) = 0.151425. 

We study the evolution of Q32 throughout run C3. We 
find that Q32 (and therefore the density perturbation 
with angular dependence given by Y32) is indeed propor- 
tional to a, at least as long as the growth of a remains 
exponential; after that, Q^2 grows more slowly than a, 
and it reaches a maximum a few rotation periods before a 
(see Fig. |l^. The phase evolution of the density pertur- 
bation is also consistent with expectations: the Fourier 
transform of Q32{t) shows a very definite peak at the 
r-mode (numerical) frequency w. 

A quantitative check shows that Eq. ( p^ ) predicts the 
observed magnitude of Q32 with an accuracy of about 
50%; this error is consistent with the next-order terms 
(^ fl'^ and a^) not included in this expression. In the 
slowly rotating calibration model CI, we find that Q32 is 
given by Eq. ( |36|) to within about 1%. 

We point out that we do not explicitly include any 
density perturbation in the initial configuration of the 
star; rather, the density perturbation is immediately gen- 
erated by the hydrodynamic evolution of the fluid as a 
consequence of the initial velocity perturbation. The evo- 



lution of the amplitude of the density perturbation am- 
plitude provides more insight into the mechanism that 
causes the r-mode to saturate: on the surface of the star, 
dp appears as four large wave crests; at a critical ampli- 
tude these crests stop growing, and within a few rotation 
periods they turn into breaking waves that damp the r- 
mode. 



H. Limits on mode—mode coupling 

In the numerical evolution (C3) nonlinear hydrody- 
namic processes do not prevent the gravitational radia- 
tion instability from driving the dimensionless amplitude 
of the r-mode to values of order unity. In particular, the 
energy of the r-mode is not channeled into other modes 
by nonlinear hydrodynamic coupling until the amplitude 
of the mode becomes quite large. It is possible however 
that the nonlinear processes that would limit the growth 
of the r-mode act only on timescales that are longer than 

(s) 

our artificially brief simulation growth time t^^, but still 
shorter than the physical r^^. 

Can our numerical simulation place any limits at all 
on the possibility of nonlinear coupling? We know that 
in our simulation the amplitude of the r-mode grows ex- 
ponentially until a w 2, so the nonlinear interaction with 
other modes must be negligible at least until that time. 
This observation allows us to set a limit on the strength 
of the nonlinear couplings between the modes; and from 
this limit we can infer a lower limit on the saturation am- 
plitude that may be achieved when the radiation-reaction 
coupling is adjusted to its physical value. Of course, the 
inference is only justified for the nonlinear interaction of 
the r-mode with other modes that are correctly modeled 
in our simulation (for instance, the finite azimuthal res- 
olution of the grid sets an upper limit on the m of the 
modes that can be resolved), and with our physical as- 
sumptions (for instance, the buoyant 5-modes of realistic 
neutron stars will not be present with our choice of the 
equation of state). 

Our argument is based on the Lagrangian description 
of the nonlinear evolution of the mode amplitude devel- 
oped by Schenk, et al. |8|. In this formalism, the modes 
interact at the lowest order by way of three-mode cou- 
plings: roughly speaking, quadratic interactions between 
pairs of modes drive the evolution of the amplitude of a 
third mode. Because at the beginning of our simulation 
all modes except the r-mode have negligible amplitude, 
we expect that the most important three-mode nonlin- 
ear term might be one that couples two r-modes to a 
third mode . Following Ref. Q we consider the coupled 
equations for the r-mode and a generic mode X obtained 
in second-order Lagrangian perturbation theory: 



R 



dc 
~dt 
dcx 
dt 



VjJXCx 



CR 

trr ^ 

iujx K*XRR 

2 ex 



2 Cfi 



'^r'^r^ 



(37) 
(38) 
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where cr and cx are the complex amphtudes (including 
phases) of the modes; and ujx are their frequencies; 
eR and ex are the nonlinear mode energies at unit am- 
plitude; and trr is the radiation-reaction e-folding time 
of the r-mode. Finally, Kxrr is the nonlinear interaction 
energy for unit amplitude modes. Schenk, et al. give 
expressions for the kxrr of coupled generic Newtonian 
modes in rotating stars. In writing Eqs. ( ^7|) and (|38| ) we 
have omitted the coupling terms proportional to k,*xxw 
which are forbidden by a z-parity selection rule ||^ : the r- 
mode has odd z parity, so it cannot couple quadratically 
to the mode X. 

From our numerical evolution C3, we know that the 
amplitude of the r-mode grows very nearly exponentially 
until a ~ 2: 



CR{t) ~ Cfl(0)e 



(39) 



(s'] 

where r^j^ is the artificially short radiation-reaction 
timescale used in our simulation. (Although it is con- 
venient to take \cr\ ~ a, our argument still applies as 
long as \cr\ is merely proportional to a.) Therefore, we 
also know that until \cr\ ~ 2, the second term on the 
right side of Eq. ( |37| ) is negligible compared to the first. 
In this case. 



^RR 



> 



(40) 



We now use Eq. (^9|) to integrate Eq. ( |38| ) and compute 
cx- 



exit) = cx{0)e 



— iujxt 



^^x >i*xRR Kit)]^ ~ [c*RiO)?e - 
2iujR + iujx + 2/- 



As) 

RR 



(41) 



Now we set cx(0) ~ and |c/j(t)| 3> |ch(0)| for the time 
late in the simulation when c/j ~ 2, and find 



\cx{t)\ 



^XBR 



ex 



^xl41l4(t)|^ 



^RR 



(42) 



where Sto = 2ujr + lox ■ We define the resonance index 
7^*) = \ujr/ujx\[{t'^^6ujY + 4]^/^, whose value is close 
to unity, 7^^^ ~ 1, when the system is near resonance, 
6uj ~ 0. We use this bound on |cx(i)| in Eq. ( pO| ) to 
obtain 



1 



> 



''XRR\ 



,(^) UxtR ^^^^ 7^'^) 

'RR ' 



(43) 



We can rewrite this inequality in terms of the r-mode 
period Pr = 2'k/lor-. 



R 



'RR 



\cR{t)\' 



ex 



ex 
efl 



(44) 



We now set |cii;(t)| = 2 (the value at which the evolution 
of the amplitude begins to show deviation from exponen- 

(s) 

tial) and Pr/t^^^ = 1/10 (the value for our simulation), 
and obtain 



'^xrr 



ex 



' ex 7^-^ 
tR 4007r2 ■ 



(45) 



Thus, our numerical evolution puts a limit on the 
strength of the coupling between the r-mode and other 
modes in the star. 

We now ask how the saturation amplitude would 
change if the radiation-reaction timescale assumed its 

iv) is) 

physical value instead of the value r^^j^ used in our 
simulation C3. The key to doing this is to realize that 
Eqs. ( |37| ) and ( |3^ ) describe the coupled mode evolution 
in the physical case if we just substitute for t'^^. 
The mode X is capable of stopping the unstable growth 
of the r-mode only when the magnitude of the second 
term on the right side of Eq. ( |37| ) becomes comparable 
to the first. Through an analysis similar to the one that 
led to Eq. (^) , it is straightforward to find the following 
condition on the saturation amplitude of the r-mode. 



Pr 



^RR. 



„sat 1 2 



7 



(p) 



^XRR 



ex 



ex 
e_R ' 



(46) 



We now use the upper limit for Ik^^^^ I from Eq. ( [iq ) from 
our numerical evolution, to obtain a lower limit for the 
amplitude c^* at which the r-mode would be saturated 
in the physical case: 



icri 



> 20- 



Pr. 



ip) 

'RR 



(47) 



Since 7^^^ > 7^*' this equation yields |c|f*| 3> 4x 10^** for 
run C3. So if the dominant mode-mode coupling is of the 
form given in Eqs. ( ^7|) and (^8|), our simulation places a 
relatively large lower limit on the r-mode saturation am- 
plitude. However, the r-mode could instead be limited by 
parametric resonance pl[ with a suitable pair of modes 
(satisfying the resonance condition ojr + LOy + ojz — 0). 
It appears that our simulation does not provide a very 
strong lower limit on the saturation amplitude that could 
be imposed by this kind of process. 



I. Dependence on the grid spacing 

We wish to confirm that our standard computational 
grid can resolve the spatial structure of the r-mode well 
enough to give reliable predictions about the saturation 
amplitude of the mode. For this purpose, we have per- 
formed a simulation (run C3*) with the same parameters 
of run C3, but on a grid with only half the spatial resolu- 
tion (i. e., 32 cells in the radial direction, and 64 cells in 
the axial and azimuthal directions). Figure compares 
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FIG. 16: Numerical evolution of the r-mode amplitude a in 
low-resolution run C3* (solid curve) and in run C3 (dashed 
curve) . 



the evolution of a in runs C3 and C3*. The two curves 
are very similar, but in run C3* saturation is reached 
a bit earher, at t/Po — 21 A, and at a somewhat lower 
amplitude a — 2.68. This may be caused by the larger 
numerical viscosity that must be present in the coarser 
grid. The evolution of the other diagnostics is also very 
similar in the two runs. 

Thus, the simulation run (C3*) suggests that the qual- 
itative results of our simulations are independent of the 
resolution adopted. The r-mode saturation amplitudes 
on the two grids agree to within about 20%. Interest- 
ingly, the extrapolation to the infinite resolution case 
suggests that the physical saturation amplitude might 
be even larger than 3.3. 



FIG. 17: Numerical evolution of the r-mode amplitude a in 
production runs C3-C5. 
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FIG. 18: Numerical evolution of the r-mode frequency to in 
production runs C3-C5. 



VI. TESTING THE SATURATION 
AMPLITUDE 

Even in the absence of a saturation mechanism due to 
mode-mode coupling as described above, it is possible 
that the saturation amplitude in our simulation might 
still depend on the strength of the radiation-reaction 
force. In our simulation we see that the r-mode grows 
until density waves on the surface of the star break and 
form shocks. It is possible that this occurs just because in 
our simulation we are pushing the fluid too hard with the 
radiation-reaction force, much harder than it would be 
appropriate in the physical case. To explore how the evo- 
lution depends on the strength of this driving force, we go 
back to the time in run C3 before any signs of nonlinear 
saturation are seen, when a = 1.8. We start a new run 
(C4) there, increasing the value of k (which determines 
the strength of the radiation- reaction force) to 5967 (1.33 
times its value in run C3). The new growth timescale 
is about 7.5Po- (Undoubtedly, a test with k <^ 4487 
would have been more compelling; but our evolutions 
are so computationally expensive that we were forced to 
increase rather than decrease the strength of the driving 



force.) 

In a separate run (C5), we test the influence of the 
history of the evolution of the r-mode on its saturation 
amplitude. Namely, we ask if an r-mode that started out 
as the linear initial data of Eq. (p3|), with a very large 
amplitude, would evolve much differently from an r-mode 
that started out small and was built up gradually to large 
amplitude by the radiation reaction force. To answer this 
question, we start with the Fast equilibrium model, and 
we add a linear r-mode velocity fleld with oq = 1.8. For 
this run we keep k — 4497. 

Figures [ij, ^ and |l^ show the evolution of the di- 
agnostic parameters a, lu and Ail for runs C3-C5. As 
expected, the r-mode does grow faster in run C4, but its 
maximum value is essentially the same (the maximum 
a = 3.338 at t = 24.12Po) as in run C3. In this run, the 
r-mode amplitude increases from a ~ 1.8 to a = 3.338 
within a time At ~ 6Po (compared to At ~ 8Po in run 
C3) as would be expected given that the driving force is 
I times that of run C3. 

In run C5, the growth of the r-mode is initially slower 
than in run C3, as the linear r-mode velocity field evolves 
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FIG. 19: Numerical evolution of the differential rotation AQ 
in production runs C3-C5. 



2.0 F 




FIG. 20: Numerical evolution of the r-mode amplitude a in 
production runs C6-C8. 



toward its correct nonlinear form. Eventually its maxi- 
mum occurs at essentially the same amplitude as before 
(a — 3.337). Figures [l7| and show that during run C5 
a and lu undergo short-period oscillations; this happens 
because the initial velocity field is only a small-amplitude 
approximation to the real m = 2 r-mode eigenfunction. 
So other spurious modes with fairly large amplitude are 
excited initially in run C5. Note that these extraneous 
modes must make nonzero contributions to J22 if they 
are to show up in our diagnostics. Here the extrane- 
ous modes cause a rapid modulation of a and lu with a 
dominant period of about O.SPq- Finally, it is interesting 
to consider the evolution of Afl (Fig. |l9|), which is very 
similar in the three runs. 

These runs provide limited evidence that the satura- 
tion amplitude of the r-mode does not depend (strongly) 
on our artificially large radiation reaction force. The non- 
linear hydrodynamical process that leads to shock forma- 
tion appears to be triggered by attaining a certain crit- 
ical amplitude of the r-mode, with little dependence on 
the strength of the radiation-reaction force. Thus if no 
mode-mode coupling occurs on timescales longer than 

(s) 

our unphysically short t)j^, then our results suggest that 
the maximum amplitude a « 3 is a reasonable guess for 
the physical case (k = 1) as well. 



VII. FREE EVOLUTION 

Stergioulas and Font 0| have also studied the nonlinear 
evolution of r-mode initial data, but using relativistic 
hydrodynamics in a fixed background geometry. In their 
evolution using this relativistic Cowling approximation, 
the gravitational interactions of the mode with itself and 
with the rest of the star are neglected. The principal 
difference between their model and ours therefore is that 
theirs has no radiation reaction, and no r-mode growth. 

Stergioulas and Font find that, for an initial r-mode 
amplitude ao = 1-0, no significant suppression of the 



mode is observed during 13 rotation periods. They define 
their mode amplitude using a post-Newtonian exp ression 
for the eigenfunction that differs from our Eq. ( |23|) except 
in the Newtonian limit. And their method of evaluating 
the mode amplitude numerically also differs from ours. 
They read the mode amplitude from the value of the 
fiuid's velocity at a single point within the star, while we 
define a in terms of integrals over the entire star. In the 
slow-rotation Newtonian limit our two definitions agree. 
Stergioulas and Font observe that the amplitude of the 
velocity oscillations (shown in Fig. 2 of Rcf. j^]) decrease 
by about 50% during the course of their simulation, an 
effect that they attribute to numerical viscosity [Q. In 
order to compare our own simulations more directly with 
theirs, we performed a series of evolutions in which we 
turned off the radiation-reaction force by setting k = 0. 

In production runs C6 and C7, we augment our rapidly 
rotating equilibrium configuration with the approximate 
r-mode velocity field of Eq. (p3|). For run C6, we choose 
the initial ao so that a [as measured by our numerical 
diagnostic, Eq. (^)] is initially 1.8: the value at which we 
start to observe deviations from exponential growth in 
run C3. For run C7, we choose ao so that the initial a 
is 1.0, in order to make a direct comparison with Ster- 
gioulas' and Font's published results. We have evolved 
these systems through respectively 11 and 7 initial rota- 
tion periods (several hydrodynamical timescales, accord- 
ing to our rough estimate of the speed of sound for the 
rapidly rotating model). 

We plot the evolution of a and lu for these simula- 
tions in Figs. and |l]. The wavy appearance of the 
curves suggests that, by using the linear eigenfunction, 
Eq. (p3[), for amplitudes of order unity, we have excited 
spurious modes in addition to the basic m = 2 r-mode. 
We have already observed this behavior in run C5. The 
rapid modulation of a and oj has a period of about O.SPq, 
and the amplitude of the modulation is smaller for run 
C7. (This is reasonable: for lower a we expect the ap- 
proximate expression, Eq. (p^), to be more accurate and 
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FIG. 21: Numerical evolution of the r-mode frequency u) in 
production runs C6 and C7. 
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start a new run (C8) using the C3 data at this time. We 
evolve these data setting k to zero in the subsequent evo- 
lution. We follow this evolution through an additional 
15.4 initial rotation periods. During this time the mode 
amplitude a is essentially constant, see Fig. except 
for a slow secular decline due to numerical viscosity at 
0.23%/revolution, and a few very small amplitude oscil- 
lations. The r-mode frequency is quite constant, and the 
phase coherence function, and the differential rotation 
Arj also remain quite small in this case (see Fig. p2l ). 
The r-mode amplitude in run C3 remains above unity 
for 14.3 rotation periods, so run C8 demonstrates that 
the LSU hydrodynamic code 10 1 used here reliably 
and stably evolves large amplitude r-modes in rapidly 
rotating stars for the duration of our simulations. 

Comparing runs C6, C7, and C8, we infer that the 
strong decrease in the amplitude observed in runs C6 
and C7 occurs as nonlinear hydrodynamics reorganizes 
the initial linear r-mode velocity field to the correct non- 
linear form for amplitudes of order unity. After the re- 
organization is complete (within a few rotation periods), 
a decreases only because of numerical viscosity. (In run 
C5, this same phenomenon caused the slower growth of 
the amplitude compared to run C3.) By contrast, the 
small decrease in run C8 appears to be caused entirely 
by numerical viscosity. 

Altogether, we find that our results are compatible 
with those of Stergioulas and Font 0| : no nonlinear sat- 
uration effect is evident in the free nonlinear evolution of 
r-modes, at least for amplitudes of order unity. 



FIG. 22: Numerical evolution of the differential rotation Af2 
in production runs C6-C8. 



SO to excite smaller amplitude spurious modes.) 

In both runs, a loses about 20% of its initial value 
during the first four rotation periods. In the next few ro- 
tation periods, however, the average value of a remains 
unchanged (although in run C6 we can see a further mod- 
ulation of the amplitude with a period of about 8Po)- 
Throughout the runs, the r-mode frequency lo oscillates 
around uj — —1. 12170, consistent with its value in run 
C3 for the same value of a (i. e., 1.44). As the run is 
started, the differential rotation All (which is zero in 
the initial, rigidly rotating star) increases almost imme- 
diately to values that are consistent with those observed 
in run C3 for the same amplitude; compare Figs. ^ and 
p^ . As a decreases, Afi decreases consistently. (In run 
C7, Ail settles to a value slightly higher than what we 
expected from its value in run C3 when a — 0.82; but we 
did not run this evolution as far as run C6, so at the end 
of our simulation the value of All might still be evolving.) 

Finally, we study the free nonlinear evolution of an r- 
mode that was grown to the amplitude a = 1. To do 
so, we go back to the time in run C3 when a = 1, and 



VIII. REPEATED SPINDOWN EPISODES? 

The first attempt to analyze the nonlinear evolution 
of r-modes by Owen, et al. Q was based on a simple 
two-parameter model consisting of a rotating star with 
angular velocity and its r-mode with amplitude a. Us- 
ing this model the mode was found to grow exponentially 
until it reached some maximum level amax, where it was 
assumed to remain saturated. Energy and angular mo- 
mentum were expected to be removed from the star by 
gravitational radiation during this saturation phase until 
the r-mode regained stability (because of increased inter- 
nal dissipation brought about by cooling or because the 
angular momentum of the star was reduced to a very low 
level). In this initial picture gravitational radiation was 
expected to spin down the star on a timcscale of about 
one year. The radiation emitted was expected to sweep 
down in frequency from | times the initial angular ve- 
locity of the star to | times its final value: ranging from 
perhaps 1 kHz initially to perhaps 100 Hz. 

Our simulations suggest a very different picture. We 
find that, once the amplitude of the r-mode reaches amax, 
it is quickly reduced by the action of the breaking waves 
and shocks, instead of remaining saturated at this value 
for a very long time. At the end of our simulation the 
star still has 60% of its initial angular momentum, and 
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FIG. 23: Numerical evolution of the r-mode amplitude a in 
the extended run C4. 



FIG. 24: Numerical evolution of the r-mode frequency oj in 
the extended run C4. 



its average angular velocity is 77.5% of JIq. Thus the 
star is left rotating relatively rapidly, leaving open the 
possibility of subsequent episodes of r-mode instability 
and spindown. 

To investigate this possibility, we extend run C4, evolv- 
ing our star for 13 more initial rotation periods after a 
has gone back to its initial value (0.1), or (equivalently) 
for nine periods after a reaches its minimum (~ 0.01). 
The evolution of the amplitude for this case is plotted 
in Fig. After t — 33Po, the fluid motion is quite 
turbulent, but we see no sign that a is starting to grow 
again. The evolution of the r-mode frequency (Fig. |2^ ) 
is also erratic, probably because here the sinusoidal ap- 
proximation begins to fail (remember that uj is approx- 
imated as —{XlJi-2)d\J'2,2\ldt^- In fact, after t = 33Po 
we have found it necessary to impose an ad hoc limit on 
the value of w; otherwise w grows to about — 17rio, and 
the radiation-reaction force (proportional to oj^) becomes 
huge, pushing the fluid to superluminal velocities. 

Nine periods should be more than enough to see a sec- 
ond r-mode growth episode, if it occurs at all. Although 
at the end of the simulation the average angular velocity 
of the star is lower than f^o, the growth timescale is de- 
termined by the r-mode frequency, which is even higher 
than at the beginning of the run. What keeps the r-mode 
then from resuming its growth? 

One hypothesis is that because of its strong differential 
rotation the post-spindown conflguration of the star is 
one which stabilizes the r-mode. The value of Af2 for the 
last few periods is plotted in Fig. ^ The increase of Afi 
observed between t = 32Po and t = 36Po is not caused by 
radiation reaction, but by a global, energy-conservative 
reorganization of the fluid. At the end of this process, the 
spatial structure of differential rotation is very different 
from what it was at Ofmax: compare Fig. |lj (t = 25.6Po 
in run C3) with Fig. |§ (t = 42Po in run C4). The latter 
plot shows a star that is rotating on cylinders (except for 
the outer layer), with ri(iz7, z) almost proportional to vj. 

Karino, et al. p3] derived linearized structure equa- 




FIG. 25: Differential rotation Af2 through the extended run 
C4. 



tions for the r-modes of differentially rotating Newtonian 
stars. When differential rotation is so strong that coro- 
tation points appear (that is, when there exists a such 
that Lu + mfl{w) = 0), the mode equations go singular. 
(The presence of a corotation point at the cylindrical 
radius vj means that the velocity pattern of the mode 
appears to stand still in the frame rotating with angular 
velocity Q^-cu).) A comparison of the differential rotation 
of Fig. with the value of to suggests the presence of 
corotation points in the final configuration of our star. 
By itself, however, the singularity of the linearized mode 
equations does not necessarily mean that r-modes are 
impossible. 

A second, probably more likely possibility is that, in 
the very noisy environment manifest in Figs. ^ and 
Figs. E4, the growing r-mode is unable to get locked in 
phase with the approximate expression for the driving 
force that we use here. The actual radiation reaction 



force [Eq. (17)] is a function of the frequency of the r- 
mode. Since we do not know exactly what this frequency 
is, we use the expression (^l]) to approximate it. This ap- 
proximation works extremely well as long as the r-mode 
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FIG. 26: Meridional structure of differential rotation at the 
end of production run C4. This contour plot shows level con- 
tours for the value of the azimuthally averaged angular veloc- 
ity Q{zz!, z)/Q, at time t — 42Po. 



makes the dominant contribution to J22 , yet, in the tur- 
bulent post-spindown environment, the r-mode no longer 
dominates the evolution of J22. Hence, our expression for 
the gravitational radiation reaction force is no longer cor- 
rect: it fails to maintain phase coherence with the r-mode 
and so prevents the growth of the mode. 

If the r-mode really does not exist in the chaotic post- 
spindown environment, then it will be necessary to wait 
for viscosity to damp differential rotation before the r- 
mode can grow again. However, viscosity might be un- 
able to do this before the star cools so much that the r- 
mode is stabilized (either because the star forms a crust 
or because viscosity itself has grown too strong). This 
possibility is worrisome, because the same environmen- 
tal conditions (strong differential rotation and general- 
ized noise) that characterize the end of run C4 are likely 
to occur in the young supernova remnants where r-modes 
are expected to arise in nature. Still, we think it more 
likely that the absence of a second growth episode in our 
simulation is the result of our expression for the radia- 
tion reaction force, which is too simple for this chaotic 
situation. 



IX. CONCLUSIONS 

We have completed a series of numerical 3D hydro- 
dynamical simulations of the nonlinear evolution of the 
GR driven instability in the r-modes of rotating neutron 
stars. We have verified that the current-quadrupole GR 



reaction force implemented in our code is accurate by re- 
producing the analytical predictions (for slowly rotating 
stars) with our full 3D numerical integration code. In 
our simulations, the amplitude of the (m = 2) r-mode 
is driven to a value of about three before nonlinear hy- 
drodynamic forces stop its growth by the formation of 
shocks and breaking surface waves. We showed that the 
value of this maximum amplitude is insensitive to the 
strength of the GR driving force by repeating the sim- 
ulation for different strengths and different initial fluid 
configurations. We also repeated our simulation using 
a coarser numerical grid to verify the robustness of our 
results (the maximum mode amplitude changes only by 
about 20% when the number of grid points is reduced by 
a factor of 8), and to show in particular that numerical 
viscosity is not playing a critical role in our simulations. 

In our simulation we have artificially increased the 
strength of the GR reaction force in order to reduce the 
problem to one that can be studied with the available 
computer resources. We have shown, however, that the 
results of our simulation can be used to infer limits on the 
real physical problem as well. We used the results of our 
simulations to derive a lower limit of a few times 10~* on 
the saturation amplitude of the r-mode in a real neutron 
star due to possible (but unseen) nonlinear mode-mode 
couplings. This lower limit applies to couplings with 
modes that are well described by our simulation: that 
is, the modes of a barotropic fluid with spatial structures 
larger than about 2% of the radius of the star. 

Recent analysis of the effects of magnetic fields [ p3[ , 
and exotic forms of bulk viscosity Q suggest that the 
r-mode instability may not play as important a role in 
astrophysical situations as was once thought. However, 
the considerable uncertainty that exists about both the 
macroscopic and microscopic states of a neutron star 
makes it impossible at the present time to conclude that 
the r-mode instability plays no astrophysical role. Thus 
it seems reasonable to us that some effort be put into 
gravitational wave searches for r-mode signals having 
forms qualitatively similar to those predicted by simu- 
lations such as this. 
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APPENDIX A: USEFUL EXPRESSIONS IN 
CYLINDRICAL COORDINATES 

In this Appendix we give explicit expressions in cylin- 
drical coordinates (tu, z, if) for a number of useful quan- 
tities used in our simulations. The components of the 
the initial r-mode velocity field used in our numerical 
evolutions are 



, 5 r!o . „ 
aoA/ Y^-^ztn sm2(p, 



5 2 • r. 



(Al) 



(A2) 



and 



5 

Q,qvj + ai)\j — — —zmcos2if. (A3) 

V lOTT it 

We refer the azimuthal component of the velocity to the 
orthonormal coordinate (p, so that v'^ and have the 
same numerical value and we can use them interchange- 
ably. 

The integrals that determine J22 and its first time- 
derivative are, 



'22 



(A4) 



(1) 



7' 

where 

Ti 
T2 



- j pe-^"^[Ti + iT2]wdzudzdip, (A5) 
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(A6) 



,9$ 



9$ 



z{v^ - v^) - wvr^Vz + 07 zm—.{A7) 

^ oz oru 



The components of the radiation-reaction force in 
cylindrical coordinates are obtained from Eq. ( p7| ) by ex- 
pressing the current multipole tensor Sjk in terms of the 
current multipole moments J2m via Eqs. (]l9|)-(^l|): 



f: 



GR 



16 IAttG 
'45' 



(A8) 



and 



45 V 5 c7 



'iVzJ22 + ^"^22^ 



(A9) 

where k = 1 in general relativity theory. The fifth 



(5) 

and sixth time derivatives of J22 are obtained as J22 — 
uj'^J^ , and = —uj^J22- 
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